#include <math.h>
#include "sleputil.h"

void slep_randn(double* x, int n, double mu, double sigma){
	static int i=0;
	const double p0 = 0.322232431088;     const double q0 = 0.099348462606;
	const double p1 = 1.0;                const double q1 = 0.588581570495;
	const double p2 = 0.342242088547;     const double q2 = 0.531103462366;
	const double p3 = 0.204231210245e-1;  const double q3 = 0.103537752850;
	const double p4 = 0.453642210148e-4;  const double q4 = 0.385607006340e-2;
	double u, t, p, q, z;
	if(i==0){
		i++;
		PlantSeeds(-1);
	}

	for(i=0; i<n ; i++){
		u   = Random();
		if (u < 0.5)
			t = sqrt(-2.0 * log(u));
		else
			t = sqrt(-2.0 * log(1.0 - u));
		p   = p0 + t * (p1 + t * (p2 + t * (p3 + t * p4)));
		q   = q0 + t * (q1 + t * (q2 + t * (q3 + t * q4)));
		if (u < 0.5)
		z = (p / q) - t;
		else
		z = t - (p / q);
		x[i] = (mu + sigma * z);
	}
}